Photon Blockade in the Ultrastrong Coupling Regime 
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We explore photon coincidence counting statistics in the ultrastrong-coupling regime where the 
atom-cavity coupling rate becomes comparable to the cavity resonance frequency. In this regime 
usual normal order correlation functions fail to describe the output photon statistics. By expressing 
the electric-field operator in the cavity-emitter dressed basis we are able to propose correlation 
functions that are valid for arbitrary degrees of light-matter interaction. Our results show that the 
standard photon blockade scenario is significantly modified for ultrastrong coupling. We observe 
parametric processes even for two-level emitters and temporal oscillations of intensity correlation 
functions at a frequency given by the ultrastrong photon emitter coupling. These effects can be 
traced back to the presence of two-photon cascade decays induced by counter-rotating interaction 
terms. 



PACS numbers: 42.50.Pq, 42.50.Ar, 85.25.-j, 03.67.Lx 

The quantum theory of photodetection and optical co- 
herence as originally formulated by Glauber pQ is central 
to all of quantum optics and has occupied a key role 
in understanding radiation-matter interactions. Pho- 
todetection is also common to quantum-state engineer- 
ing [3j and quantum information protocols [3]. One of 
the most prominent effects that shows nonclassical be- 
havior of a quantum emitter is a phenomenon known as 
photon blockade @H6], where a coherent excitation of a 
cavity coupled to a highly nonlinear quantum system, 
such as a single atom, quantum dot or superconduct- 
ing qubit, generates an output train of single photons. 
This photon statistics for the cavity output can be inves- 
tigated by measuring the intensity correlation function 
g^(r), which demonstrates the nonclassical character of 
the transmitted field. 

Recently a new regime of cavity quantum electrody- 
namics (cavity QED) has been reached experimentally 
where the strength of the interaction between an emitter 
and the cavity photons becomes comparable to the tran- 
sition frequency of the emitter or the resonance frequency 
of the cavity mode [THT2] . In this so called ultrastrong 
coupling regime pj [HI I13H16] the routinely invoked ro- 
tating wave approximation is no longer applicable. As 
a consequence the number of excitations in the cavity- 
emitter system is no longer conserved, even in the ab- 
sence of drives and dissipation. 

In this letter we explore the photon statistics of the 
output fields of a driven cavity QED system where a two 
level system (TLS) interacts with one cavity mode in 
the ultrastrong coupling regime. Specifically we show 
that the standard photon blockade effect does no longer 
persist since ultrastrong emitter-photon couplings enable 
processes where a single photon that enters the cavity 
can generate two or more output photons. Whereas such 
parametric processes can of course occur for emitters 



with three or more levels, we here show their existence 
for a two level system. We find that the two Rabi-splitted 
resonances, which are often termed upper and lower po- 
lariton, exhibit very different photon statistics. While 
excitation of the lower energy peak provides output light 
which is both subpoissonian and antibunched, as in the 
standard regime, excitation of the higher energy peak 
results in the emission of superpoissionan and bunched 
light. The calculated resonance fluorescence spectrum 
shows that this result originates from the activation of 
second order nonlinear optical effects which are absent 
in usual two level atoms. Such a process can only result 
from counter-rotating terms in the atom-cavity interac- 
tion Hamiltonian. 

As a second main result we find that the time depen- 
dence of photon-photon correlations g^ (r) differs signif- 
icantly from the standard photon blockade scenario. For 
standard photon blockade, one finds 5 (2) (fJ) < 1 and os- 
cillations with the Rabi frequency of the laser drive for 
t > (provided ft exceeds the loss rates). The latter 
emerge since a time ~ is needed to load the next 
photon into the nonlinear cavity once a photon has left. 
For the ultrastrong coupling regime, we find oscillations 
at a frequency of the order of the emitter-photon cou- 
pling g ^> fl, i.e. much faster oscillations. These emerge 
since upper and lower polariton excitations can not be 
generated independently. 

Our findings can be measured in present day experi- 
ments [7l412j . A particularly well suited technology for 
such an experiment are superconducting circuits [T71 [TB] 
which have recently emerged as an exquisite platform for 
microwave on-chip quantum-optics experiments. Even 
though single photon detectors in the microwave regime 
are under current development [El [20], Glauber's first 
and second order correlation functions have been mea- 
sured [3TJ [52] using quadrature amplitude detectors and 
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standard photon blockade at microwave frequencies [33] 
has been observed. 

Input-output relations - Applying Glauber's idea of 
photodetection, we here introduce a full quantum theory 
to study the photon blockade in the ultrastrong coupling 
regime. This requires a proper generalization of input- 
output theory 24J, since the standard relations would 
for example predict an output photon flux that is pro- 
portional to the average number of cavity photons, i.e. 
(do Ut (t)a ou t(t)) oc (a^(t)a(t)) for vacuum input. Hence an 
incautious application of this standard procedure to the 
ultrastrong-coupling regime would predict an unphysi- 
cal stream of output photons for a system in its ground 
state which contains a finite number of photons due to 
the counter-rotating terms in the Hamiltonian. To rem- 
edy this problem, a model which explicitly takes into 
account the colored nature of the dissipation bath has 
been proposed [32] US] • Yet such a procedure is numeri- 
cally quite demanding and would require the derivation 
(and also the time-integration) of four-time correlation 
functions for the calculation of g^ 2 \r). 

Here we propose a different route going back to 
Glauber's formulation of photodetection [1 J. By ex- 
pressing the cavity electric-field operator in the atom- 
cavity dressed basis we derive correlation functions for 
the output fields which are valid for arbitrary degrees of 
light-matter interaction in the cavity. The probability 
per unit time that a photon be absorbed by an ideal 
detector is proportional to (E~(t)E + (t)) were E {t) 
are the positive and negative frequency components of 
the electric field operator of the output field [TJ [37J 
and higher order photon correlation funtions e.g. read 
{E-{t)E-{t')E+{t')E+{t)). In circuit QED the same 
quantities are measured via output voltages which are 
proportional to the electric fields. 

To derive the appropriate input-output relations '2A\, 
we assume a cavity that is coupled to a one-dimensional 
output waveguide via an interaction between the cavity 
field X and the momentum quadratures II W of the waveg- 
uide field outside the cavity, Hi = e c J duXH^ with 

n w = -iyj [a u (t) - al(t)] , where e c is a coupling 
parameter and e Q is a parameter describing the dielectric 
properties of the output waveguide, v is the phase ve- 
locity and a ul (a[ J ) annihilation(creation) operators of the 
fields outside the cavity. This form of the coupling ap- 
plies to optical implementations where the electric fields 
in- and outside the cavity couple as well as to circuit 
QED setups where voltages (and hence electric fields) of 
the cavity and output line couple via a capacitance. For 
optical realizations, e G is thus the permittivity and for cir- 
cuit QED setups the capacitance per unit length of the 
output line. We define output (input) field operators, 

«out(i„)(i) - -j= J dutfje-W-Qao®, (1) 



where t — t\ — >• +oo for the output field and i — 
to —> — oo for the input field. With this defini- 
tion, h(al ut (t)a ou t(t)) is proportional to the measured 
(E~(t)E + (t)) and describes an energy flux associated to 
the output light. The definition of output fields as used 
in many textbooks, c.f. [35], is recovered if all frequencies 
of the field are very close to a "carrier" frequency cU and 
one may approximate y/uj rs \/W in the integral kernel 
which makes the observed energy fluxes proportional to 
photon number fluxes. Following the standard procedure 
we obtain the input-output relation, 

flout (t) = flin(i) ~ i /o % X+, (2) 

where we have applied a rotating wave approximation 
to Hj since uj >• \e c /{\/&K I Kejv)\. One thus needs to 
find the positive frequency component of X according 
to its actual dynamical behavior, c.f. [29], to compute 
the proper output fields. We do this by expressing X 
in the atom-cavity dressed basis. Importantly, in the ul- 
trastrong coupling regime, the positive frequency compo- 
nent of X is not proportional to the photon annihilation 
operator a. We now turn to introduce our model. 

Model - We consider a coherently driven cavity QED 
system with the most general linear coupling between 
a single cavity mode and a two level system (TLS). Its 
Hamiltonian reads (we set % = 1), 

H s = Ho + H dlivc with (3) 
Hq = u}$a} a + uj K <7 + a~ + g{a + a}) [cos(#)ct z — sin(#)cr x ] 

Here, H is the energy of the TLS described by the Pauli 
operators er x ,y, z (c* = (o- K ±ia y )/2) and the cavity mode 
with annihilation (creation) operators a(cv). g describes 
the strength and the mixing angle 9 the further prop- 
erties of the coupling between the TLS and the cavity 
mode. An additional coupling term cx (a + a))o~ y would 
not change the physics as it could be compensated by a 
rotation around the z-axis. -ffdrivc = ^lcos(u>^t)(a + (v) 
describes the coherent driving of the cavity mode with 
frequency Wd and drive amplitude f2. The mixing an- 
gle 9 plays a crucial role since both, the spectrum of H n 
and the output photon statistics strongly depend on its 
value. In particular, for 9 = it/ 2 the Hamiltonian Hq 
conserves the parity of the number of excitations in the 
system, whereas this is no longer the case for 9 ^ tt/2 
[3"Ull32| . The latter enables transitions that are forbid- 
den in the usual Jaynes-Cummings (JC) model, c.f. Fig. 
[l] and leads to novel effects in the photon statistics. 
Note that for weaker couplings (g -C cj ,w x ) the term 
g(a + <v) cos(9)a z can be neglected in a rotating wave 
approximation and parity is always conserved. 

For describing a realistic system, dissipation induced 
by its coupling to the environment needs to be consid- 
ered. Yet, owing to the very high ratio g/uio, a standard 
quantum optical master equations fails as it would for ex- 
ample predict that even zero temperature environments 
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FIG. 1: (color online) Energy spectrum of Ho as function of 
the coupling strength for 9 — n/2 (left panel) and 9 — 0.93 
(right panel). In both panels lj x — loo- The vertical line marks 
the coupling strength (<? = 0.2 wo) that is used in the sequel. 
For 9 = 7r/2 the level structure is analogous to that of the JC 
model. The arrows indicate possible transitions of radiative 
decay. 



could drive the system out of its ground state. A viable 
description of the system's coupling to its environment 
requires a perturbative expansion in the system-bath cou- 
pling strength. To accurately perform this expansion we 
write the Hamiltonian in a basis formed by eigenstates 
\j) of Hq, denote the respective energy eigenvalues by 
Tiojj, i.e. H \j) = hujj\j}, and derive Redfield equations 
[55] to describe the dissipative processes [51]. We choose 
the labeling of the states \j) such that u>k > loj for k > j 
and focus on a single-mode cavity and a T — tem- 
perature environment. Generalizations to a multi-mode 
cavity and T ^ environments are straightforward. As- 
suming a weak laser drive, VL <C g,u>o,u; x , and treating CI 
perturbatively to leading order, see supplementary infor- 
mation, we thus arrive at the master equation [39 , 



p(t) = i[p(t),H s ] + C a p(t) + C x p{t). 



(4) 



Here C a and C x are Liouvillian superoperators de- 
scribing the losses of the system where C c p(t) = 
T, j , k>j n h W){k\W) for c = a,a~ and V[G\p = 
\{20pO^ - pO^O - O^Op). Standard dissipators are 
recovered in the limit g — > 0. The relaxation coefficients 
Tl k = 2ird c (A kj )a 2 (A kj )\Cj k \ 2 depend on the spectral 
density of the baths d c (A k j) and the system-bath cou- 
pling strength a c (A k j) at the respective transition fre- 
quency Afcj — uo k — ujj as well as on the transition co- 
efficients Cjk — —i(j\(c — c>)\k) (c = a, a~). These re- 
laxation coefficients can be interpreted as the full width 
at half maximum of each \k) —¥ \j) transition. As we 
consider a cavity that couples to the momentum quadra- 
tures of fields in one-dimensional output waveguides, we 
assume the spectral densities d c (A k j) to be constant and 
a^.(Akj) oc Afcj. Hence the relaxation coefficients reduce 
to T{ k — j c \C^ k \ 2 , where j c are the standard damp- 
ing rates of a weak coupling scenario. In equation Q 
contributions of dephasing noise were omitted as these 
only lead to very minor quantitative modifications of our 



FIG. 2: (color online) Emission characteristics of our system, 
a) (/ 2 '(0) as function of the frequency, wa, of the coherent 
drive for 9 — 0.93 (red continuous line) and 9 — n/2 (blue 
dashed line), b) Resonance fluorescence spectra calculated 
for a driving field tuned with the transition |0) — >■ i.e. the 
lower polariton, (green dashed line) and with the transition 
0) — > 1 2), i.e. the upper polariton, (orange dash-dotted line) 
for 9 = 0.93. The remaining parameters are lo x = ojq, 7 a = 
Jx = 10~ 2 loq, = 10 -4 too, g = 0.2 loo for all plots. 



findings, see supplementary information. 

Results - According to the input-output relation ^ 
the photon statistics of the output light is, for input fields 
in vacuum, equal to, 



(X-(t)X-(t + T)X+(t + T)X+(t)) 

*^oo (X-(t)X+{t)) 2 



where X = —iX a (a — a r ). We thus compute g( 2 '(r) from 
the stationary state solution of Eq. ([ij. To separate X 
in its positive and negative frequency components, X + 
and X~ , we expand it in terms of the energy eigenstates 
\j) and find X+ = -iJ2j,k>j & kj X jk \j)(k\, where X jk = 
(j\X\k) and X- = (X+)t. Note that X+\0) = 0, for the 
system ground state |0) in contrast to a\0) ^ 0. 

Fig. [2j a), shows g^(r — 0) for lo x = uj , 7 a = 
7^ = 10~ 2 w , fi = 1CT 4 lj , g = 0.2 ljq. The blue- 
dashed line represents </ 2 )(0) for the parity conserving 
model (9 = 7r/2), where a suppression of the proba- 
bility to measure two photons per count is evident for 
the two polaritonic peaks. This has the standard ex- 
planation that only one photon at a time can be ab- 
sorbed on the respective transitions due to the nonlin- 
earity of the system. Instead, the strong bunching be- 
tween both polaritonic peaks (ujd ~ loq) is due to the 
simultaneous excitation of both polaritons. By varying 
8 we observe a new and anomalous effect. For 9 = 0.93 
the expected anti-bunching for the higher frequency dip 
(cj(j s» 1.18u>o) disappears and turns into a pronounced 
bunching (# (2) (0) ~ 5.6), c.f. FigJ|a). Wc find a similar 
behavior for lower values of 9. 

An explanation of this exotic behavior can be found 
from the dynamics of the decay and excitation mecha- 
nisms, where it is crucial which transitions between en- 
ergy levels exhibit radiative decay. To probe the latter 
we calculated the incoherent part of the scattered light, 
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i.e. the resonance fluorescence spectrum which reads 
S{uj) oc lim 23? / (X~{t)X + {t + T ))e^ T dT (6) 

t->oo J Q 

for input fields in vacuum. Here 5i denotes the real part. 
For a weak excitation density S(u>) is proportional to the 
emission. Fig. [2] b) shows the resonance fluorescence 
spectrum for a weak driving field resonant with the tran- 
sitions |0) -> |1) (blue dashed line) and |0) -> |2) (red 
continuous line). The first case exhibits a single peak of 
typical Lorenzian shape reminiscent of a weakly driven 
TLS. Here, only level |1) is efficiently excited due to 
strong non-linearity and hence emission predominantly 
occurs on the transition |1) — > |0). Instead, when the 
system is resonantly excited on the |0) — > |2) transition, 
the resonance fluorescence displays a triplet structure. 
These peaks arise from the \2) —> |0), |2) — ^ |1) and 
|1) — > |0) transitions. This triplet is a peculiarity of the 
interaction terms proportional to o~ z in H , for which the 
matrix elements Xf k for all involved levels j, k = 0, 1, 2 
are non-zero and thus enable parametric multi-photon 
processes with probabilities proportional to |A^| 2 . Here 
we create an excitation in level 2 by resonantly driving 
the |0) — > 1 2) transition and the system can decay via the 
cascade \2) — s- |1) — > |0) which gives rise to the bunch- 
ing effect observed in <?( 2 )(0), c.f. Fig. [2] a). Whenever 
9 = it/2, the standard dipolar selection rules as in the JC 
model are recovered and only transitions between states 
that differ by one excitation number are allowed. Hence, 
there is no radiative decay on the \2) — > |1) transition and 
g( 2 ) (0) shows the characteristic antibunching even if we 
drive the system resonantly with |0) — > \2) (blue dashed 
line in Fig. |ja)). 
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FIG. 3: (color online) g ' (t) for a driving frequency resonant 
to the energy of the first (blue dashed line) and on the second 
polariton (red continuous line). o; x = wo, j a = "fx — 10~ 2 cjo, 
Q = 10~ 4 ujo, g = 0.2 ujo and 8 = 0.93 for both cases. 

Moreover, driving the |0) — > |2) transition triggers os- 
cillations of </ 2 ^(t) between bunching and antibunching 
values at a frequency A 2 i as shown in Fig. [3j Such oscil- 
lations also appear in a driven Jaynes-Cummings system 
albeit at much smaller amplitude but are absent for a 
single driven two level atom. Their appearance can be 
explained by the presence of different types of excitations. 
For a two level atom there is only one type of excitations 



with one resonance frequency in the system. Hence all 
dynamics in (r) is due to loading and emission pro- 
cesses only. For a driven Jaynes-Cummings system, there 
are two types of excitations, the lower polaritons and the 
upper polaritons which differ in frequency by A21 ~ 2g 
and one sees oscillations of this frequency in g^{r). Yet 
if one drives one of the Rabi peaks, only one polariton 
species is predominantly excited and the amplitude of 
the oscillations is very small. In contrast, if one drives 
the system of Eq. ^ in the ultrastrong coupling regime 
on the transition |0) — > |2), one always significantly ex- 
cites two excitations species whose frequencies differ by 
A21 which in turn causes oscillations of large amplitude 
in <t 2 '(t). Oscillations of the same frequency also ap- 
pear in g^ 2 \r) for w d = uj 01 , c.f. Fig. [3] (blue dashed 
line), but are much smaller in amplitude as the level |3) 
is barely excited. The appearance of these pronounced 
oscillations can thus also be traced back to parametric 
processes enabled by the ultrastrong coupling strength. 

Circuit QED implementation - We finally discuss an 
experimental setup with superconducting circuits that is 
ideally suited for observing our findings [35] . The Hamil- 
tonian of Eq. ([3| is realized in a coplanar waveguide 
that is galvanically coupled to a flux qubit [8]. In this 
implementation, the mixing angle 9 is related to the bias 
flux through the qubit loop. For ultrastrong couplings 
a description with a TLS and a single cavity mode is a 
valid approximation provided the flux threaded through 
the qubit loop is chosen such that cos 9 < 0.48, see sup- 
plementary information. 

Conclusions - We have explored the output photon 
statistics of a driven cavity QED system that operates 
in the ultrastrong coupling regime. We find that the 
prominent photon blockade does not persist in its known 
form for such strong couplings since parametric processes 
emerge. For efficient coupling between cavity and output 
modes, the latter could be exploited for building high 
yield parametric down-converters |36| . Generalizations of 
our study to frequency resolved detection [37| and multi- 
cavity devices [38] would form interesting perspectives 
for future research. 

This work is part of the Emmy Noether project HA 
5593/1-1 and was supported by the CRC 631, both 
funded by the German Research Foundation, DFG. 
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SUPPLEMENTARY INFORMATION 
Perturbative treatment of the driving term 

A derivation of a quantum master equation leads to 
the equation [24], 

p(t) = -i[H s (t),p s (t)] (7) 

- f dsTr E {[HT,[Hr(s-t),pE®p{t)]\}, 

J to 

where H r (s -t) = U(t - s)HjW(t - s) with, 

U(t - s) = e - iH «( { - s ) T+e- l ti~ s dTHs(r) (8) 

Here Hs(t) is the time dependent Hamiltonian of our 
system that includes the time dependent drive term, c.f. 
equation (3) of the main text, 7+ is the appropriate time 
ordering operator and He is the Hamiltonian of the en- 
vironment. In turn, p(t) is the density matrix of the 
system and pe the state of the environment which is the 
vacuum in our case. Our perturbative treatment of the 
drive consists in approximating, 



U(t-s) 



e -iH E (t~s) e -iHo{t-s) 



(9) 



where Hq is the system Hamiltonian without the drive 
term as in the main text. For, f2 -C g,tdo,uj x this is a 
good approximation since the right hand side of equation 
([7]) is already 2nd order in the weak system environment 
coupling. Our approach is thus a second order expansion 
in both, the system environment coupling and the drive 
amplitude f2. 

We furthermore apply a rotating wave approximation 
which discards all processes where system and environ- 
ment respectively system and drive would both gain 
(loose) energy at the same time. 

Due to the time dependence of H$ (t) , even the asymp- 
totic state, lim t _ i . 00 p(t), has small residual oscillations. 
For calculating steady state expectation values we aver- 
age in time over several such residual oscillations. This 
closely corresponds to the time integration of the out- 
put signal that is usually applied to record experimental 
data. 



Circuit QED implementation 

Experimental setups with superconducting circuits are 
ideally suited for observing our findings. The Hamilto- 
nian of Eq. (3) in the main text is realized in a trans- 
mission line resonator that is galvanically coupled to a 
flux qubit [HJ. In this implementation, the mixing an- 
gle 8 is related to the bias flux through the qubit loop. 
With respect to the experimental feasibility for measur- 
ing the effects we predict, it is important to verify the 
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TABLE I: 





2.25 GHz/630 nA 


27T / 27T / 2TT 


2.782 GHz/5.357 GHz/7.777 GHz 


31 1 92 1 33 
2tt 1 In 1 2-rr 


314 MHz/636 MHz/568 MHz 



validity of the description with a two level system (TLS) 
and a single cavity mode. The nonlinearity of typically 
employed flux qubits [8] greatly exceeds the qubit transi- 
tion frequency so that a TLS description is very accurate. 
As we consider the lowest frequency mode of the cavity, 
couplings to higher modes can be neglected provided that 
processes that de-excite the qubit and destroy a photon 
in the lowest frequency mode to create a photon in the 
second harmonic mode are ineffective. We now turn to 
estimate when this is the case. 

The Hamiltonian for the setup comprising all resonator 
modes reads, 



n 

+ X! ( a i + a n )(cos8(T z ~ sm6(j x ) 



(10) 



with, cos 9 = 



21 „ 



^/A2 + (2I p 8<i>y 



and sin# 



y/Ai + (2I p 84>) 2 



where uj q and ui n are the qubit and resonator frequencies 
respectively, g n is the qubit resonator-mode coupling and 
5(f), A and I p are the flux threaded through the loop, 
the sweet spot frequency 0J q (S(j) = 0) and the persistent 
circulating current of the flux qubit. Numerical values for 
all these parameters extracted from [5] are tabularized in 
table [I] The value of the qubit transition frequency as 
used in the main text is the value of uj q for the chosen 
flux bias 54>o, he. wx = w g (<$</>o). 

Our goal is to find parameter ranges where the em- 
bedded flux qubit despite the utrastrong coupling regime 
only couples to the fundamental, i.e. the lowest frequency 
resonator mode. To this end we consider the eigenstates 
of a noninteracting system of flux qubit and resonator 
which are direct products of the eigenstates of the qubit 
and the resonator, e.g. ni^n^) and have eigenener- 
gies E = ±hu) q +J2i hwi(ni + ^). Here, frequency shifts of 
the qubit and resonator due to their mutual presence are 
already taken into account, c.f. [35] . Experiments show 
that it suffices to take only the three lowest eigenmodes 
of the resonator into account Bj . The qubit-resonator in- 
teraction couples these eigenstates. For the single-mode 
description of the cavity to be valid, couplings between 
states Itlj^iCbOa) and states |tl)0in,2n3) with n 2 > 1 or 
n-3 > 1 need to be negligible. They can indeed be ne- 
glected provided the coupling is small compared to the 
difference of their energies. The eigenenergies of the non- 
interacting system are plotted as a function of the flux 
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FIG. 4: Eigenenergies of the uncoupled resonator and flux 
qubit system as a function of the flux threaded through the 
qubit loop. The energies of states that belong to the Hilbert 
space of our description (TLS and lowest frequency resonator 
mode) are drawn in blue whereas states that do not belong 
to this Hilbert space are drawn in red. Anticrossings due to 
interaction are labeled by the order of perturbation theory 
that lifts the degeneracy. The crosshatched area marks the 
allowed range of S</>- 



threaded through the qubit loop in figure [4] We focus on 
the area around the sweet spot (8(f> = 0) and thus need 
to keep Sip small enough so that the energy difference be- 
tween the state |f I1O2O3) and the state \\. O1I2O3) stays 
large compared to their mutual coupling. The energies of 
states |i W1O2O3) with m > 2 have been omitted in figure 
[4] since they are either well separated from \\. O1I2O3) or 
only cross \X O1I2O3) for much larger values of 8(f). 

On resonance between |f li0 2 3 ) and U-O1I2O3) we 
find for their coupling, 



. gi<?2 cos 9 sin 9 



UJ2 



(11) 



in second order perturbation theory in 171 and 172 ■ Despite 
the ultrastrong coupling, perturbation theory yields a re- 
liable estimate for the values of g/1^0 = 0.2 as used in the 
main text. 

We need to make sure that the energy difference be- 
tween the two states always exceeds this coupling energy 
J. This restricts the allowed range of the flux threaded 
through the qubit loop and in turn cos 9 is bound from 
above by cos 9 max « 0.48 



Effect of qubit dephasing 

In typical circuit QED implementations of the ultra- 
strong coupling regime, qubit dephasing can not be ne- 
glected. To analyse its impact on our findings, we present 
here results for a system where dephasing takes place at 
the same rate as damping, i.e. 7deph = 7- 

We model the dephasing with a Liouvillian superop- 
erator £ deph p(t) = r^ cph 2?[| j) (j|]p(i), where r^ cph = 
7dc P h0Vz|j) and V[0]p = \{WpO^ - pO^O - O^Op), 
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FIG. 5: (color online) Comparison between <?^ 2 '(0) as cal- 
culated without dephasing noise (red continuous line) and 
g^ 2 \0) as calculated with dephasing noise (green dash-dotted 
line) for 9 — 0.93. oj x = cjo, ^ = 10~ 4 ljo, g — 0.2 ujq, 
7 a = 10 -2 ujo and 7^ = 10~ 2 ujq (red continuous line) re- 
spectively 7de P h = "fx = 0.5 x 10~ 2 wo (green dash-dotted 
line). 



c.f. [31]. Figure [H] shows a comparison between </ 2 )(0) as 
calculated without dephasing noise (red continuous line) 
and g^(0) as calculated with dephasing noise (green 
dash-dotted line) for 6 = 0.93. uj x = uj 0} £1 = 10 -4 uj , 
g = 0.2 cjq, 7 a = 10~ 2 ujq and 7^ = 10~ 2 uiq (red contin- 
uous line) respectively 7d e ph = lx = 0.5 x 10~ 2 uj (green 
dash-dotted line). Hence the red continuous line in fig- 
ure [5] is exactly the same as the red continuous line in 
figure 2a of the main text. In particular for the drive fre- 
quencies of interest close to the polariton resonances, the 
differences between both results are vanishingly small. 



